/********************************************************************************
* Install HonestDiD
local github https://raw.githubusercontent.com
net install honestdid, from(`github'/mcaceresb/stata-honestdid/main) replace
honestdid _plugin_check

********************************************************************************
* Install 
* 	- coefplot (for plots, obviously)
*	- ftools (unclear if necessary; includes required functions?)
*	- reghdfe (unclear if necessary; commands also appear to work with xtreg)
*	- plot scheme (aesthetic preference)

local github https://raw.githubusercontent.com
ssc install coefplot,      replace
ssc install ftools,        replace
ssc install reghdfe,       replace
net install scheme-modern, replace from(`github'/mdroste/stata-scheme-modern/master)
net install parallel, from(https://raw.github.com/gvegayon/parallel/stable) replace
mata mata mlib index
set scheme modern

********************************************************************************/

* NLS 
use "${data}\NLS_hh_data.dta", clear

* Set up treatment and timing variables
gen mdate = ym(calendar_year,calendar_month)
local t0 = ym(2020,2)
gen int mig_covid = cond(migr_any, mdate, `t0')
format mdate mig_covid %tm

* Grab dates for graph lines and axis labels
* Min and max dates
qui: summ mdate
local tmin = r(min)
local tmax = r(max)

* COVID-19 date
local cdate = 721.5

* Season dates
* Boro harvest
local yrmin = 2018
local yrmax = 2020

local ub = 0.4
local lb = -0.12

forvalues yr = `yrmin' / `yrmax' {
	local bmin`yr' = ym(`yr',4) /* Start date Apr 1 */
	local bmax`yr' = ym(`yr',5) /* End date May 1 */
	local boro_ub`yr' scatteri `ub' `bmin`yr'' `ub' `bmax`yr'', recast(area) color(navy%15) lw(none)
	local boro_lb`yr' scatteri `lb' `bmin`yr'' `lb' `bmax`yr'', recast(area) color(navy%15) lw(none)
	if `yr' == `yrmin' {
		local boro (`boro_ub`yr'') (`boro_lb`yr'')
	}
	else {
		local boro `boro' (`boro_ub`yr'') (`boro_lb`yr'')
	}
}
di "`boro'"

* Aman harvest
forvalues yr = `yrmin' / 2019 {
	local amin`yr' = ym(`yr',11) + 15 / 30  /* Start date Nov 15 */
	local amax`yr' = ym(`yr',12) + 10 / 31  /* End date Dec 10 */
	local aman_ub`yr' scatteri `ub' `amin`yr'' `ub' `amax`yr'', recast(area) color(navy%15) lw(none)
	local aman_lb`yr' scatteri `lb' `amin`yr'' `lb' `amax`yr'', recast(area) color(navy%15) lw(none)
	if `yr' == `yrmin' {
		local aman (`aman_ub`yr'') (`aman_lb`yr'')
	}
	else {
		local aman `aman' (`aman_ub`yr'') (`aman_lb`yr'')
	}
}
di "`aman'"

* Lean season
forvalues yr = `yrmin' / 2019 {
	local lmin`yr' = ym(`yr',9)  /* Start date Sep 1 */
	local lmax`yr' = ym(`yr',11) + 10 / 30  /* End date Nov 10 */
	local lean_ub`yr' scatteri `ub' `lmin`yr'' `ub' `lmax`yr'', recast(area) color(maroon%15) lw(none)
	local lean_lb`yr' scatteri `lb' `lmin`yr'' `lb' `lmax`yr'', recast(area) color(maroon%15) lw(none)
	if `yr' == `yrmin' {
		local lean (`lean_ub`yr'') (`lean_lb`yr'')
	}
	else {
		local lean `lean' (`lean_ub`yr'') (`lean_lb`yr'')
	}
}
di "`lean'"

* Baseline DiD
reghdfe hungry b`t0'.mig_covid [aw=wgt_c], a(hhid mdate) cluster(hhid) noconstant

* Create matrix for plotting coefficients on continuous x-axis
preserve
keep mig_covid
sort mig_covid
duplicates drop
mkmat mig_covid, mat(months_t)
mat months = months_t'
mat list months
restore

* Coef plot and export
local plotopts ///
	title("") graphregion(color(white)) plotregion(margin(sides)) ///
	mcolor(gs0) ciopts(recast(rcap) color(gs0)) 
graph set window fontface "CMU Serif"
coefplot, vertical  omitted baselevels ///
	at(matrix(months)) xline(`cdate', lcolor(gs0) lpattern(shortdash)) text(0.41 `cdate' "COVID-19") ///
	xlabel(`tmin'(1)`tmax', format(%tmMon_yy) angle(45) nogrid) ///
	 ytitle("Estimate and 95% Conf. Int.") ///
	 ylabel(-.1(.1).4, glcolor(gs12%25) glsty(solid)) ///
	 addplot(`boro' `aman' `lean', below) ///
	`plotopts'
graph export "$figures\Figure_S1.png", hei(2400) wid(3200) replace

********************************************************************************
* HonestDiD with relative magnitude restrictions

* Baseline: all pre periods, mean of post period coefficients
matrix l_vec_base = 0.5 \ 0.5

honestdid, pre(1/18) post(20/21) mvec(0.5(0.5)2) coefplot l_vec(l_vec_base) ///
	parallel(4) alpha(0.05) ///
	ytitle("Confidence Interval for Treatment Effect (95 percent)") ///
	xtitle("Mbar") ///
	`plotopts' 
graph export "${figures}\Figure_S2.png", hei(2400) wid(3200) replace

* Baseline: "p-value" for M = 1
honestdid, pre(1/18) post(20/21) mvec(1) l_vec(l_vec_base) ///
	parallel(4) alpha(0.124)

/* OTHER SPECS 
* All pre periods, first post period coefficient
matrix l_vec_0 = 1 \ 0
honestdid, pre(1/18) post(20/21) mvec(0.5(0.5)2) coefplot l_vec(l_vec_0) ///
	parallel(4) `plotopts' alpha(0.05)
graph export "${output}\nls_honestdid_mr_all_0.png", hei(2400) wid(3200) replace

* All pre periods, second post period coefficient
matrix l_vec_1 = 0 \ 1
honestdid, pre(1/18) post(20/21) mvec(0.5(0.5)2) coefplot l_vec(l_vec_1) ///
	parallel(4) `plotopts' alpha(0.05)
graph export "${output}\nls_honestdid_mr_all_1.png", hei(2400) wid(3200) replace
*/

* Re-estimate using separate month and year FEs
reghdfe hungry b`t0'.mig_covid [aw=wgt_c], a(hhid calendar_month calendar_year) cluster(hhid) noconstant

coefplot, vertical  omitted baselevels ///
	at(matrix(months)) xline(`cdate', lcolor(gs0) lpattern(shortdash)) text(0.41 `cdate' "COVID-19") ///
	xlabel(`tmin'(1)`tmax', format(%tmMon_yy) angle(45) nogrid) ///
	 ytitle("Estimate and 95% Conf. Int.") ///
	 ylabel(-.1(.1).4, glcolor(gs12%25) glsty(solid)) ///
	 addplot(`boro' `aman' `lean', below) ///
	`plotopts'
graph export "${figures}\Figure_S3.png", hei(2400) wid(3200) replace

honestdid, pre(1/18) post(20/21) mvec(0.5(0.5)2) coefplot l_vec(l_vec_base) ///
	parallel(4) alpha(0.1) ///
	ytitle("Confidence Interval for Treatment Effect (90 percent)") ///
	xtitle("Mbar") ///
	`plotopts' 
graph export "${figures}\Figure_S4.png", hei(2400) wid(3200) replace

/* OTHER SPECS
* Re-estimate regression using only matching months
reghdfe hungry b`t0'.mig_covid if inlist(month,2,3,4) [aw=wgt_c], a(hhid mdate) cluster(hhid) noconstant

* Matching months, mean of post period coefficients
honestdid, pre(1/6) post(8/9) mvec(0.5(0.5)2) coefplot l_vec(l_vec_base) ///
	parallel(4) `plotopts' alpha(0.05)
graph export "${output}\nls_honestdid_mr_mon_mean.png", hei(2400) wid(3200) replace

* Matching months, first post period coefficient
honestdid, pre(1/6) post(8/9) mvec(0.5(0.5)2) coefplot l_vec(l_vec_0) ///
	parallel(4) `plotopts' alpha(0.05)
graph export "${output}\nls_honestdid_mr_mon_0.png", hei(2400) wid(3200) replace

* Matching months, second post period coefficient
honestdid, pre(1/6) post(8/9) mvec(0.5(0.5)2) coefplot l_vec(l_vec_1) ///
	parallel(4) `plotopts' alpha(0.05)
graph export "${output}\nls_honestdid_mr_mon_1.png", hei(2400) wid(3200) replace
*/
********************************************************************************
* NPL
use "${data}\NPL_hh_data.dta", clear
replace round_date = td(01jun2020) if round==7	// to look pretty

* Set up treatment and timing variables
drop if round <= 5
local t0 = date("2019-07-16","YMD")
gen int high_covid = cond(highRemit, round_date, `t0')
format round_date high_covid %td

* Grab dates for graph lines and axis labels
* Min and max dates
qui: summ round_date
local tmin = r(min)
local tmax = r(max)

* COVID-19 date
local cdate = date("2020-03-05","YMD")

* Create matrix for plotting coefficients on continuous x-axis
preserve
keep round_date
sort round_date
duplicates drop
local rounds = _N
mkmat round_date, mat(dates_t)
mat dates = dates_t'
mat list dates
restore

* Appendix regression to match
areg FoodIns treat covid i.month, a(hhno) vce(cluster hhno)
reghdfe FoodIns treat covid, a(month hhno) cluster(hhno)
reghdfe FoodIns treat, a(round_date hhno) cluster(hhno)

* Season dates
* Rice harvest
local yrmin = 2018
local yrmax = 2020

local ub = 2
local lb = -1.15

forvalues yr = 2018 / 2019 {
	local rmin`yr' = date("`yr'-10-11","YMD") /* Start date Oct 11 */
	local rmax`yr' = date("`yr'-11-15","YMD") /* End date Nov 15 */
	local rice_ub`yr' scatteri `ub' `rmin`yr'' `ub' `rmax`yr'', recast(area) color(navy%15) lw(none)
	local rice_lb`yr' scatteri `lb' `rmin`yr'' `lb' `rmax`yr'', recast(area) color(navy%15) lw(none)
	if `yr' == `yrmin' {
		local rice (`rice_ub`yr'') (`rice_lb`yr'')
	}
	else {
		local rice `rice' (`rice_ub`yr'') (`rice_lb`yr'')
	}
}
di "`rice'"

* Wheat harvest
forvalues yr = 2019 / `yrmax' {
	local wmin`yr' = date("`yr'-04-08","YMD") /* Start date Apr 8 */
	local wmax`yr' = date("`yr'-05-15","YMD") /* End date May 15 */
	local wheat_ub`yr' scatteri `ub' `wmin`yr'' `ub' `wmax`yr'', recast(area) color(navy%15) lw(none)
	local wheat_lb`yr' scatteri `lb' `wmin`yr'' `lb' `wmax`yr'', recast(area) color(navy%15) lw(none)
	if `yr' == `yrmin' {
		local wheat (`wheat_ub`yr'') (`wheat_lb`yr'')
	}
	else {
		local wheat `wheat' (`wheat_ub`yr'') (`wheat_lb`yr'')
	}
}
di "`wheat'"

* Lean season
forvalues yr = `yrmin' / 2019 {
	local lmin`yr' = date("`yr'-08-02","YMD")  /* Start date Aug 2 */
	local lmax`yr' = date("`yr'-10-11","YMD") /* End date Oct 11 */
	local lean_ub`yr' scatteri `ub' `lmin`yr'' `ub' `lmax`yr'', recast(area) color(maroon%15) lw(none)
	local lean_lb`yr' scatteri `lb' `lmin`yr'' `lb' `lmax`yr'', recast(area) color(maroon%15) lw(none)
	if `yr' == `yrmin' {
		local lean (`lean_ub`yr'') (`lean_lb`yr'')
	}
	else {
		local lean `lean' (`lean_ub`yr'') (`lean_lb`yr'')
	}
}
di "`lean'"

* X-axis labels
local xlabels = ""
local mon_min = mofd(`tmin')
local mon_max = mofd(`tmax')
forvalues tick = `mon_min' / `mon_max' {
	local dt = dofm(`tick')
	local xlabels = "`xlabels' `dt'"
}
di "xlabels is `xlabels'"

* Coef plot and export
reghdfe FoodIns b`t0'.high_covid, a(hhno round_date) cluster(hhno) noconstant

local plotopts ///
	title("") graphregion(color(white)) plotregion(margin(sides)) ///
	mcolor(gs0) ciopts(recast(rcap) color(gs0))

coefplot, vertical omitted baselevels ///
	at(matrix(dates)) xline(`cdate', lcolor(gs0) lpattern(shortdash)) text(2.06 `cdate' "COVID-19") ///
	xlabel(`xlabels', format(%tdMon_yy) angle(45) nogrid) ///
	 ytitle("Estimate and 95% Conf. Int.") ///
	 ylabel(-1(.5)2, glcolor(gs12%25) glsty(solid)) ///
	 `plotopts' ///
	 addplot(`rice' `wheat' `lean', below)
graph export "${figures}\Figure_S5.png", hei(2400) wid(3200) replace

* Parameter to average sensitivity across post-COVID periods
matrix l_vec_base = 0.5 \ 0.5

* HonestDiD with relative magnitude restrictions
honestdid, pre(1/11) post(13/14) mvec(0.5(0.5)2) coefplot l_vec(l_vec_base) ///
	ytitle("Confidence Interval for Treatment Effect (95 percent)") ///
	xtitle("Mbar") ///
	parallel(4) `plotopts' alpha(0.05)
graph export "${figures}\Figure_S6.png", hei(2400) wid(3200) replace

* Baseline: "p-value" for M = 1
honestdid, pre(1/11) post(13/14) mvec(1) l_vec(l_vec_base) ///
	alpha(0.477)

gen mon = month(round_date)
gen year = year(round_date)
reghdfe FoodIns b`t0'.high_covid, a(hhno mon year) cluster(hhno) noconstant

coefplot, vertical omitted baselevels ///
	at(matrix(dates)) xline(`cdate', lcolor(gs0) lpattern(shortdash)) text(2.06 `cdate' "COVID-19") ///
	xlabel(`xlabels', format(%tdMon_yy) angle(45) nogrid) ///
	 ytitle("Estimate and 95% Conf. Int.") ///
	 ylabel(-1(.5)2, glcolor(gs12%25) glsty(solid)) ///
	 `plotopts' ///
	 addplot(`rice' `wheat' `lean', below)
graph export "${figures}\FIgure_S7.png", hei(2400) wid(3200) replace

honestdid, pre(1/11) post(13/14) mvec(0.5(0.5)2) coefplot l_vec(l_vec_base) ///
	ytitle("Confidence Interval for Treatment Effect (90 percent)") ///
	xtitle("Mbar") ///
	parallel(4) `plotopts' alpha(0.1)
graph export "${figures}\Figure_S8.png", hei(2400) wid(3200) replace
